/* benchmark guru interfaces */
#ifndef QMCPLUSPLUS_P3DFFT_BENCH_H
#define QMCPLUSPLUS_P3DFFT_BENCH_H
#include <complex>
#include <vector>
#include <fft/p3dfft.h>

namespace APPNAMESPACE {

template<typename T>
struct P3DFFT_benchmark
{
  typedef  T               real_type;
  typedef  std::complex<T> complex_type;

  enum {D=3};

  ///default constructor
  P3DFFT_benchmark():ng_global_r(0),ng_local_r(0),ng_local_k(0) {}

  ///constructor with the size information
  P3DFFT_benchmark(int *dist_dim,int nx, int ny, int nz)
    :ng_global_r(0),ng_local_r(0),ng_local_k(0) 
  {
    resize(dist_dim,nx,ny,nz); 
  }

  ///destructor
  ~P3DFFT_benchmark()
  {
    p3dfft_clean();
  }

  ///return the size of the global data
  inline size_t get_global_size_r() const { return  ng_global_r; }
  ///return the size of the local data
  inline size_t get_local_size_r() const { return  ng_local_r; }
  ///return the size of the local data
  inline size_t get_local_size_k() const { return  ng_local_k; }

  /** transform forward
   * @param in starting address of the input data
   * @param out starting address of the output data after the transformation
   */
  void transform_forward(T* in, T* out)
  {
    p3dfft_ftran_r2c(in,out);
  }

  /** transform forward
   * @param in starting address of the input data
   * @param out starting address of the output data after the transformation
   */
  void transform_backward(T* in, T* out)
  {
    p3dfft_btran_c2r(in,out);
  }

  void resize(int *dist_dim,int nx, int ny, int nz)
  {
    //already initialized with the same dimension
    if(ng_global_r && nx==ngrid[0]&&ny==ngrid[1]&&nz==ngrid[2]) 
      return;
    //cleanup to create new p3dfft
    if(ng_global_r)  p3dfft_clean();

    ngrid[0]=nx;
    ngrid[1]=ny;
    ngrid[2]=nz;

    //p3dfft_setup(dist_dim,nx,ny,nz,1);
    p3dfft_setup(dist_dim,nz,ny,nx,1);
    get_dims(istart,iend,isize,1);
    get_dims(fstart,fend,fsize,2);

    ng_local_r=isize[0]*isize[1]*isize[2];
    ng_local_k=fsize[0]*fsize[1]*fsize[2];
    ng_global_r=ngrid[0]*ngrid[1]*ngrid[2];
  }

  template<typename IT>
  void initialize(IT first)
  {
    const T pi = atan(1.0)*4.0;
    const T twopi = 2.0*pi;
    T facx=twopi/static_cast<T>(ngrid[0]);
    T facy=twopi/static_cast<T>(ngrid[1]);
    T facz=twopi/static_cast<T>(ngrid[2]);
    std::vector<T> sinx(isize[0]);
    std::vector<T> siny(isize[1]);
    std::vector<T> sinz(isize[2]);
    for(int i=0,ii=istart[0]-1; i<isize[0]; ++i,++ii) sinx[i]=sin(ii*facx);
    for(int i=0,ii=istart[1]-1; i<isize[1]; ++i,++ii) siny[i]=sin(ii*facy);
    for(int i=0,ii=istart[2]-1; i<isize[2]; ++i,++ii) sinz[i]=sin(ii*facz);

    for(int i=0; i<isize[0]; ++i)
      for(int j=0; j<isize[1]; ++j)
      {
        T sinxy=sinx[i]*siny[j];
        for(int k=0; k<isize[2];++k) *first++ = sinz[k]*sinxy;
      }
  }

  ///global grid size  in the real space
  size_t ng_global_r;
  ///local grid size  in the real space
  size_t ng_local_r;
  ///local grid size  in the  spectral space
  size_t ng_local_k;
  int ngrid[D];
  int istart[D],isize[D],iend[D];
  int fstart[D],fsize[D],fend[D];
};
}
#endif
